ANOMALOUS SHOCK DISPLACEMENT PROBABILITIES FOR A 
PERTURBED SCALAR CONSERVATION LAW 
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Abstract. We consider an one-dimensional conservation law with random space-time forcing 
and calculate using large deviations the exponentially small probabilities of anomalous shock profile 
displacements. Under suitable hypotheses on the spatial support and structure of random forces, 
we analyze the scaling behavior of the rate function, which is the exponential decay rate of the 
displacement probabilities. For small displacements we show that the rate function is bounded 
above and below by the square of the displacement divided by time. For large displacements the 
corresponding bounds for the rate function are proportional to the displacement. We calculate 
numerically the rate function under different conditions and show that the theoretical analysis of 
scaling behavior is confirmed. We also apply a large-deviation-based importance sampling Monte 
Carlo strategy to estimate the displacement probabilities. We use a biased distribution centered on 
the forcing that gives the most probable transition path for the anomalous shock profile, which is 
the minimizer of the rate function. The numerical simulations indicate that this strategy is much 
more effective and robust than basic Monte Carlo. 
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1. Introduction. It is well known that nonlinear waves are not very sensitive 
to perturbations in initial conditions or ambient medium inhomogeneities. This is in 
contrast to linear waves in random media where even weak inhomogeneities can affect 
significantly wave propagation over long times and distances. It is natural, therefore 
to consider perturbations of shock profiles of randomly forced conservation laws as 
rare events and use large deviations theory. The purpose of this paper is to calcu- 
late probabilities of anomalous shock profile displacements for randomly perturbed 
one-dimensional conservation laws. We analyze the rate function that characterizes 
the exponential decay of displacement probabilities and show that under suitable hy- 
potheses on the random forcing they have scaling behavior relative to the size of the 
displacement and the time interval on which it occurs. 

The theory of large deviations for conservation laws with random forcing is an 
extension of the Freidlin-Wentzell theory of large deviations (TUl [7] to partial differen- 
tial equations. This has been carried out extensively [UEIIIH] and we use this theory 
here. We are interested in a more detailed analysis of the exponential probabilities of 
anomalous shock profile displacements, which leads to an analysis of the rate function 
associated with large deviations for this particular class of rare events. We derive up- 
per and lower bounds for the exponential decay rate of the small probabilities using 
suitable test functions for the variational problem involving the rate function. This 
is the main result of this paper. The applications we have in mind come from uncer- 
tainty quantification |13| in connection with simplified models of fiow and combustion 
in a scramjet. Numerical calculations of the exponential decay rates from variational 
principles associated with the rate function have been carried out before [HI US] ■ We 
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carry out such numerical calculations here and confirm the scaling behavior of bounds 
obtained theoretically. We use a gradient descent method to do the optimization nu- 
merically and we note that its convergence is quite robust even though the functional 
under consideration is not known to be convex. This robustness suggests the Monte 
Carlo simulations with importance sampling using a change of measure based on the 
minimizer of the discrete rate function is likely to effective. Our simulations show 
that indeed such importance sampling Monte Carlo performs much better than the 
basic Monte Carlo method. 

The paper is organized as follows. In Section [2] we formulate the one-dimensional 
conservation law problem. In Section |3] we state the large deviation principle and 
identify the rate function which we will use. In Section |4] we give a simple, explicitly 
computable case of shock profile displacement probabilities that can be used to com- 
pare with the results of the large deviation theory. Section[5]contains the main results 
of the paper, which as the upper and lower bounds for the exponential decay rate of 
the displacement probabilities, under different conditions on the random forcing. We 
identify scaling behavior of these probabilities, relative to size of the displacement 
and the time interval of interest. In Section [6] we introduce a discrete form for the 
conservation law and the associated large deviations, and calculate numerically the 
displacement probabilities from the discrete variational principle. In Section [7] we 
implement importance sampling Monte Carlo based on the minimizer of the discrete 
rate function and we compare it with the basic Monte Carlo method. We end with a 
brief section summarizing the paper and our conclusions. 

2. The perturbed conservation law. We consider the scalar viscous conser- 
vation law: 

ut + {F{u)), = {Du,),, X e M, t e [0, cx)), (2.1) 
u[t = 0, x) = uo(x), a; G M. (2.2) 

Here F e C^(M) and the initial condition satisfies Ui^{x) — >■ m± as a; — > ztcx) where 
u_ < u+. We are interested in traveling wave solutions of the form U{x — jt), where 
the profile U satisfies U{x) — >■ u± as a; — >■ ±oo and the wave speed 7 is given by 
the Rankine-Hugoniot condition 

,= ^(-+)-^("^). (2.3) 

A traveling wave U with wave speed 7 exists provided the following conditions are 
fulfilled: 

F{u)- F{u-) > -f{u-u^), Vue(M+,M_), (2.4) 
F'{u+) < 7 < F\u^). (2.5) 

The first condition is the Oleinik entropy condition [M] and the second one is the Lax 
entropy condition 16J. Under these conditions the traveling wave profile exists, it is 
the solution of the ordinary differential equation 

= ^{F{U) ~ F(u_) - 7([/ - u^)), U{x) u^, 

and it is orbitally stable, which means that perturbations of the profile decay in time, 
and thus initial conditions near the traveling wave profile converge to it. Note that 
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a physically admissible viscous profile must have the stability property; otherwise, 
it would not be observable. As noted in [15 , the motivating idea behind the orbital 
stability result is that in the stabilizing process, information is transferred from spatial 
decay of the profile U at infinity to temporal decay of the perturbation. 

The purpose of this paper is to address another type of stability, that is, the 
stability with respect to external noise. We consider the perturbed scalar viscous 
conservation law with additive noise: 

ul + {F{u^))^^{Dul)^+sW{t,x), xeR, te[0,oo), (2.6) 
u%t = 0,x) = U{x), xeR. (2.7) 

Here e is a small parameter, W(t, x) is a zero-mean random process (described below), 
and the dot stands for the time derivative. We would like to address the stability of the 
traveling wave U{x — jt) driven by the noise sW . Motivated by an application mod- 
eling combustion in a scramjet [131 125] . we have in mind a specific rare event, which 
is an exceptional or anomalous shift of the position of the traveling wave compared 
to the unperturbed motion with the constant velocity 7. 

We consider mild solutions which satisfy (denoting u'^(t) — {u'^ (t, x))j:eR) 

u'{t) = S{t)U- f S{t-s)N{u'{s))ds + e f S{t-s)dW{s), (2.8) 
Jo Jo 

where S{t) is the heat semi-group with kernel 

1 / (x — y)"^ \ 

and N{u){x) = {F{u{x)))x- The main result about the heat kernel Chap. XVI, 
Sec. 3] is as follows. 

Lemma 2.1. /// = (/(i))te[o,T] e L^{[0,T],H-^{R)), then the function J^Sit- 
s)f{s)ds is in L'^{[0,Tl H^{R)) nC{[0,T], L^(R)). 

A white noise or cylindrical Wiener process B(t, x) in the Hilbert space L^(M) is 
such that for any complete orthonormal system {fn{x))n>i of L^(M), there exists a 
sequence of independent Brownian motions (/3„(i))„>i such that 

00 

B(i,x) = ^/3„(i)/„(x). (2.9) 

Tl=l 

B{t, x) can be seen as the (formal) spatial derivative of the Brownian sheet on 
[0, 00) X M, which means that it is the Gaussian process with mean zero and co- 



variance E[_B(t, x)B{s, y)] = s At 6{x — y). Note that the sum ( 2.9 1 does not converge 
in LP' but in any Hilbert space H such that the embedding from Ir(M) to H is Hilbert- 
Schmidt. Therefore, the image of the process i? by a linear mapping on L^(IR) is a 
well-defined process in the Sobolev space H^iR) (for fc = 0, we take the convention 
= L"^) when the mapping is Hilbert-Schmidt from L'^{R) into iJ'=(M). Then, if the 
kernel $ is Hilbert-Schmidt in the sense that 

00 k 
n=l j=0 
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then the random process W = ^B, 

W{t,x)^ I ^x,x')B{t,x')dx' 



is well-defined in H'^ (M) . It is a zero-mean Gaussian process with covariance function 
given by 

E[W{t, x)W{t', x')] ^tAt' C{x, x'), 

C{x,x') = <P<p'^{x,x')^ I <^{x,x")^{x' ,x")dx\ 



where stands for the adjoint of $. 

As an example, we may think at $(a;, x') = ^q{x)^i{x — x') with $o G H^{^) 
and $1 S i?'^(R). In this case $o characterizes the spatial support of the additive 
noise and $1 characterizes its local correlation fun ction. In the limit $0(2^) = 1 and 



$i(a; — x') = 5{x — x'), the process W{t, x) in (2.61 is a space-time white noise. 

By adapting the technique used in [5, Chap. XVI, Sec. 3] we obtain the following 
lemma. 

Lemma 2.2. // $ is Hilbert- Schmidt from L"^ into , then Z{t) — S{t — 
s)dW{s) belongs to L^([0,T], H^(R)) r]C{[0,T], L^{R)) almost surely. If ^ is Hilhert- 
Schmidt from into , then Z{t) belongs to L'^{[Q,T],H'^{M)) f^C{\Q,T],H^{M)) 
almost surely. 

Proof. See Appendix |Xj □ 

3. Large deviation principle. In this section we state a large deviation princi- 
ple (LDP) for the solution {u'^ {t,x))t£[o^T\,x£VL of the randomly perturbed scalar con- 



servation law (2.6). It generalizes the classical Freidlin-Wentzell principle for finite- 
dimensional diffusions. Throughout this paper, we assume that the flux _F is a 
function with bounded first and second derivatives. 

Assumption 3.1. In this paper, the flux F is a function and there exists 
Cp < 00 such that \\F'\\i^^i^-g^ and ||-F"||l°°(r) o'^e bounded by Cp. 

Remark. Assumption |3 1 1 is merely a technical assumption to simplify the proofs 
of Proposition |3.2| and |3.3| and obviously it violates the convexity of fluxes in conser- 
vation laws. From the physical point of view, for a general flux F , we can choose a 
very large constant M and let Fm{u) = F{u) for |u| < M and saturate Fm{u) for 
|u| > M. If [— Af, M] can cover the range of interest of u, then Fm{u) — F(u). With- 
out a proof, we poin t ou t that this physical argument can be proven mathematically: 



if uf, u% and in (2.6 ) are continuous on [0, T] x M and W is bounded on [0, T] x M, 
then the parabolic maximum principle implies that is bounded on [0,T] x M and 
thus we can find M < 00. The technical difficulty is that W is not differentiable in 
time so we can not apply the parabolic maximum principle directly. However, we 



can consider a modified by replacing by a smooth W in (2.6) and show that 

We first describe the functional space to which the solution of the perturbed 
conservation law belongs. 

Proposition 3.2. If ^ is Hilbert- Schmidt from to H^, then for any e > 



there is a unique solution to (2.8) in the space £^ almost surely, where 



= {{uit,x))te[o.T]..m ■ {u{t,x) - U{x))te[o.T]..eR & Ci[0,T],H\R))}. 
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Proof. See Appendix [BTTJD 

The rate function of the LDP is defined in terms of the mild sohition of the control 
problem 

Ut + {F{u)), = {Du,), + <t>h{t,x), xGR, ie[0,T], (3.1) 
u{t = 0,x) ^U{x), x£R, (3.2) 

where h G L^([0, T], L^(M)). The solution u to this problem is denoted by H[h] and 
H is a mapping from ^^([0, T], L^(M)) to £^ provided $ is Hilbert-Schmidt from 
to H^. The following proposition is an extension of the LDP proved in li. 

Proposition 3.3. //$ is Hilbert-Schmidt from to H^, then the solutions 
satisfy a large deviation principle in with the good rate function 
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T 



^(^) = . ^"£..,9 / \Ht,-)\\hdt, (3.3) 

h,u='H[h] I Jq 

with the convention inf = oo. 
Proof See Appendix [B?2l □ 

In fact, I{u) < oo if and only if u is in the range of H. The LDP means that, for 
any A (Z £^ , we have 

-J {A) < liminf e2logP(u'' e A) < limsupe^ logP(M'' e A) < - J (A), (3.4) 

with 

J{A) ^ inf I{u). (3.5) 

ueA 

Note that the interior A and closure A are taken in the topology associated with 

Although the convexity of the rate function / is unknown, it is possible to show 
that / satisfies the maximum principle: if the set of the rare event A does not contain 
the exact solution U{x — jt), then mfueA attains its minimum at the boundary 
of A. 

Proposition 3.4. IfU{x — jt) ^ A andu e A, then there exists a sequence {u"} 
in £^ such that I{u") < I{u) and u" u in £^ as n ^ oo. As a consequence, any 
u G A can not be a local minimizer of I . 

Proof. See Appendix |B.3| □ 

4. Wave displacement from an elementary point of view. In this paper we 
are interested in estimating the probability of large deviations from the deterministic 
path U{x — jt). In this section, we first study in a very elementary way how the center 
of the solution it^ at time T can deviate from its unperturbed value. The center of a 
function u G £^ is defined as 

JZ[u{t,x)-U{x)]dx 
C[u]{t) = , (4.1) 

provided that the integral is well-defined. 

Proposition 4.1. // the covariance function C is in L^(M x M), then the center 
of is well-defined for any time t G [0,T] almost surely and it is given by 

r., . iW{t,x)dx , , 

C\u^]it) ^ -it ^ ' . (4.2 



6 



It is a Gaussian process with mean jt and covariance 

Cov(CK](t),C[zi^](t')) - ,2 nC{x,x')dxdx' ^ ^ 

Proof. See Appendix [C] □ 

In the absence of noise, the center of the solution increases linearly as ^t. In the 
presence of noise we can characterize the probability of the rare event 

B = {ue£\C[u]{T)>jT + xo}, (4.4) 

where xq € [0, oo). 

Proposition 4.2. // the covariance function C is in L^iM. x M), then 

The approximate equality means that 

lim £2logP(.^ e B) = - ,/.°^"7~!|/, 
£->-o 21 jj C[x,x')axdx' 

This proposition is a direct corollary of Proposition |4.1| 

If = then U{x ~ jt) e B and F^u" e B) = 1/2. If xq > then the set B 
is indeed exceptional in that it corresponds to the event in which the center of the 
profile is anomalously ahead of its expected position. Note that the scaling Xg/T in 



(4.5) corresponds to that of the exit problem of a Brownian particle. 

The LDP for stated in Proposition |3.3| is not used here for two reasons: 

1. It does not give a good result because the interior of i? in f ^ is empty (since 
we can construct a sequence of functions bounded in that blows up in 

2. The distribution of the center C[u^]{T) is here explicitly known for any e. 
This is fortunate and it is not always true. When the rare event is more 
complex, than the LDP for is useful as we will show in the next section. 

5. Wave displacement from the large deviations point of vievif. 

5.1. The framework and result. In this section, we use the large deviation 
principle to compute the probability of the rare event that the perturbed traveling 
wave is the same profile but with the displacement xq at time T: 

A = {ue£^ such that u{0,x) = U{x), u{T,x) = U{x - jT-xq)}. (5.1) 

The value of the rate function heavily depends on the kernel $ and formally the 
simplest kernel we can have is the identity operator. Of course, the identity operator 
is not Hilbert-Schmidt, but with the given A, we can construct a Hilbert-Schmidt $ 
such that $ is approximately the identity in the region of interest. 

Assumption 5.1. 



1. We assume that the center of transition ofU{x) in (5.1) is and Xq + ^T > 0. 

2. The kernel has the following form: 

where a, Lq and 1^ are positive constants. 



7 



3. (j>o and (f>i are n functions so that their Fourier transforms are well- 
defined in M. and $ is Hilbert- Schmidt from to (see Lemma 5.4-)- 4>o 
and (f>i are normalized so that '/>o(0) — (f)i{0) — 1. 

4- (t^o G '■^ positive valued and (pi is nonzero. In addition, l/0o(a;) and 
have at most polynomial growth at ±cx). 

5. 4'oix) is increasing on x G (— oo, —1), decreasing onx £ (1, oo) and identically 
equal to 1 on x & [—1,1]. By this setting, the support of the noise kernel $ is 
roughly [-Lo,io]- 

The following theorem is the main result of this section, and the proof will be 
given in the next two subsections. 



for all satisfying Assumpti 



5.1 



Theorem 5.2. Let A be defi ned by (5.1). There exists a constant Cq such that 
-Ja „„i.- A „„„.™„-(.- — cr 1 — + xq + C'o, the quantity J^{A) of 



the large deviation problem (3.1) generated by <i>^^ has the following asymptotic (in 
the sense that Ic ^ 0) scales: 



JiA) = inf /(u) '^=° 



e{xl), 
eixim, 

0(kol), 

e(kol/r), 

e(i/T), 



for \xo\ .small and any fixed T, 
for \xo\ .small and T small, 
for \xo\ large and any fixed T, 
for \xq\ large and T small, 
for T small and any fixed xq. 



(5.3) 



Here J^{A) Q[xq) and similar expressions mean that there are constants Ci and 
C2 with the units such that Cix^ and C2XQ are dimensionless and CiXq < J^{A) < 
C2XI as Ic — > 0. 

Proof See Section [5]2] and [5^ □ 

Remark. Cq = maxjCi,, C,„}, where Cy and Cw will be determined in Lemma 
5.6| and |5.10[ respectively. Lq — 7T + xq + Cq means that the range of noise covers 
the region of interest of A so that the rare event A is possible and at the same time 
the range is not too large to cause the waste of energy; the small Ic means that the 
noise in this region is weakly correlated and can be treated as the white noise. 



We note that there are two occurrences of T in J^{A): in (3.3), we integrate the 



\\h{t, •)I!l2 from to T, and the rare event A is T-dependent. These two occurrences of 
T would make the T-dependence of J^{A) nontrivial; if T is small, the T-dependence 
in ( |5.1| is negligible and we can have the scale in T. Therefore, if the traveling wave 
is stationary (7 = 0), then A has no T-dependence so the sharper bounds can be 
obtained. 

Corollary 5.3. If j — 0, then (5.3|) can be refined: 



inf I(u) 



'Q{xl/T), 

0(ko|), 
0(kol/r), 



for \xq\ small, 

for \xq\ large and any fixed T , 
for \xq\ large and T small, 
for any fixed xq . 



(5.4) 



Proof. See Section 5.2 and 5.3 □ 

Remark. The case that both |xo| and T are large is not clear even if 7 = 0. 
This is because J{A) = e(|a;o|/T) for l^ol > T, and J{A) = Q{xl/T) for |a;o| < T- 
when both xq and T are large, J {A) is the mixture of Q{\xq\/T) and Q{x'^/T) so it 
is difficult to estimate the bounds. 
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In the rest of this subsection we discuss two relevan t iss ues about this framework. 
We first show that the given kernel in Assumption 5.1 is indeed Hilbert-Schmidt 
from to . 

Lemma 5.4. //(/)o a"'^ 0i ^"^6 both in n H'' for k > I, then the kernel of 
the form (5.2) is Hilbert-Schmidt from to . 



Proof. See Appendix D.l □ 



The other issue is that A has no interior and therefore J^{A) = oo, which gives a 
trivial lower bound in the large deviation framework; to avoid this triviality, a rigorous 
way is to consider instead the (closed) rare event Ag with a small 5 > Q: 

As = {ue £^ such that u(0, •) = [/, ||u(r, •) - [/(• - 7T - xo)||hi < S}. (5.5) 

By Proposition |3.3[ we have 



- J{Ai) < liminf e^logP(u'' e Ag) < limsupe^ logP(M'' e Ag) < -J{As). (5.6) 

In addition, the following lemma shows that J^{As) converges to J^{A) as (5 — >■ 0. 

Lemma 5.5. By definition J^{As) is a decreasing function with 6 and bounded 
from above by JT{A). In addition, 

hm J(A,) = J(A). 

(5->-0 



Proof See Appendix |D.2| □ 

By Lemma 5.5 and the fact that J^{As) < J^{A), 



'J (A) < --J{As) < liminf£^logP(w" £ Ag 



<lmisiipe-'logP{u'' e As) <-J{As) -JiA). 

e-s-O 



Namely, for e and S small, we have 



(u"^ e As) ^ exp 



tJ{A) 



and it thus suffices to consider the bounds for J^{A). 



Remark. We can compare the results stated in Theorem 5.2 valid under As 



sumption |5.H and the ones stated in Proposition |L5] First we can check that A C B. 
Second we can anticipate that in fact B is not much larger than A because it seems 
reasonable that the most likely paths that achieve C[w](T) > 'jT + xo should be of the 
stable form U{x — 7T — xq) at time t — T. This conjecture can indeed be verified as 
we now show. Note that 



C{x, x')dxdx' — 



^{x,x')dx dx' . 



If the kernel is of the form (|5.2|), then we find 



C(x, x')dxdx' ~ 



27r 



\M^O\'\MOM- 
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By assuming Ic ^ Lq (and (/)i(0) — 1) as in Assumption 5.1 then this simphfies to 



C{x^x')dxdx' ~ <T^ Lq / (f)Q{x)'^dx 



Therefore, from Theorem [5.2| (in which Lq — jT + xq + Cq) , we get 
C(x, x')dxdx' 



which gives with (4.5) 



e(xo), 

9(1), 



{u' e B) ^exp -^J{B) 



for Xq large and T small or fixed, 
for Xq and T small or fixed. 



with 



'Q{xl), 
&{xl/T), 
J{B)= {Q{\xq\), 

e(|xo|/T), 
e(l/T), 



for |a;o| small and any fixed T, 
for |a;o| small and T small, 
for |a;o| large and any fixed T, 
for |a;o| large and T small, 
for T small and any fixed ccq. 



as in (5.3). We therefore recover the same asymptotic scales as in Theorem 5.2 



5.2. Proof of Theorem 



5.2 



for small |a;o|. Here we prove the bounds for 
J{A) in (5.3) and (5.4) when |a;o| is small and the bound 0(l/r) for fixed xq. We 
first consider the upper bounds. By definition, J{A) < I{u) for any u ^ A\ therefore 
the idea to obtain the upper bound is to find good test functions u. 
For \xq\ small, we use the linear shifted profile as the test function: 

v{t,x) - U{x~{t/T){-iT + xq)). 

If the kernel w cr^^ in the region of interest, then 



J(A).in^/(.)</(«)<^ 



\\vt+F{v)^-Dv.. 



xx\\L' 



:dt. 



(5.7) 



We show (5.7) rigorously by the following lemmas: 

Lemma 5.6. Given the traveling wave U{x) in (5.1), there exists Cy > 1 such 
that, uniformly in xq and T and for all t G [0, T], 



\\[4>q\-/Lq) - IpA- - it/T)i^T + xq))\\' 



< 1. 



(5.; 



where Lq = jT + xq + d, . 

Proof. See Appendix |D.3| □ 



Lemma 5.7. Let vjt.x) = U{x - (t/T)hT + xq)). Given Lq ^ xq + + Cy in 



Lemma 



then 



5.6: satisfying Assumption 5.1 and h£^ defined by 



Vt + F{v), 



(5.9) 



\htitr)\\hdt 



h^o 1 



\\{vt + F{v)^ - Dv,^)/M-/Lo)\\l2dt. (5.10) 
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Proof. See Appendix |D.4| □ 

Now we are ready to prove (5.7). By Lemma 5.7 we have 



I{v) <\£ \\h^£,{t, -nhdt '^4° ^ 1^ IIK + F{v), - Dv,,)/M-/Lo)\\hdt. 
Then the rest is to compute \\{vt + F{v)x — Dvx-j^ I I Lof^jjidt. Note that 



With Lo = 7T + xq + Cy given in Lemma |5.6[ we have 

||C/,(- - {t/T){jT + xo))/M-/Lo)\\h < 2\\U,,{-~{t/T){jT + xomh+2, 
and therefore 

T 

\\{vt + F{v)x - Dvxx)/M-/Lo)\\l2dt 







x, 

y2 



2 fT 








^,2 pT 



\Ux{- ~ it/T){jT + XQ))/M-/Lo)\\l2dt 
<^J^ {2\m--it/T)ijT + xomh+2)dt^f{2m\l.+2). 
Therefore, we have the asymptotic upper bounds for J (A): 

J{A)<I{v)^< —^{2\\U.A\h+2). (5.11) 

Now we find the lower bomids for J {A). Let us first denote by 1 the function 
identicahy equal to 1, assume that ("I>^ )-^l e i^, and find the general form of the 
lower bound. 

Proposition 5.8. Given (<i>^^)^l e L^, for any u e A, 



I{u) > ^r-i|i(<i>^;)^l||-? (^J [U{x-x„)-Uix)]dxJ 



(5.12) 



Proof. See Appendix |D.5| □ 

Then we show that ($^^)^1 is indeed in i^. 



Lemma 5.9. Let satisfy Assumption 5.1 Then ) l{x) is in L (R), and 

ll('J^Lo)^l|li^-^'^'^o||<^o|li. as 1, ^0. 
Proof. See Appendix |D.6| □ 

From Lemma |5.6[ Lq = + + C„, where C„ is uniform in xq and T. Then 



for Ic small the general lower bound (5.121 becomes: 

J{A) > lT-'m'lfl\\-J (^J [Uix-xo)-Uix)]dxJ (5.13) 
;\\MjT-\xo + + Cy)-' (^j [U{x-XQ)-U{x)]dx^ . 
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For \xo\ small, U{x - xq) = U{x) - xqUx{x) + 0{xl). Then 



1 

2^ 



[-xaU^{x) + 0{xl)]dx] (5.14) 



\\UrjT-\xo +jT + C,)-'[xl{u+ - u^f + 0{xl)]. 



To summarize this subsection, given A in (5.11, we choose the kernel $^ with 
Lq = 7T + xq + Ci, , which is the support of the noise covering the region of interest 
of A, and then by (5.11) and (5.14), J {A) has the following asymptotic (in the sense 
that Ic — > 0) bounds: 



'Q{xl), 



^e(l/T), 
for nonzero 7, and for 7 0. 

.0 \e{xllT) 



for |a:o| small and any fixed T, 
for |a;o| small and T small, 
for T small and any fixed Xq, 



J{A) 

5.3. Proof of Theorem 



[e(i/r). 



for \xq\ small, 
for any fixed Xq. 



5.2 



for large |a;o|. Now we consider the bounds for 
J{A) in (5.3) and (5.4) when |a;o| is large. We use the test function 



w{t, x) = {l- t/T)U{x) + {t/T)U{x -jT- xo), 
the linear interpolation of the two profiles, and show that 

1 



JiA) ^ MJ{u) < liw) < 



\wt + F{w)^ - Dwxx\\^2dt. 



(5.15) 



(5.16) 



We also show (5.16) by two similar technical lemmas. 



Lemma 5.10. Given the traveling wave U{x) in (5.1) and w{t^x) in (5.15), there 
exists Cw > 1 such that, uniformly in xq and T and for all t G [0, T], 



||[0o H-/io) - l][U{- ~lT- xo) - U]\\l. < 1, 



(5.17) 



where Lq = 7T + xq + Cw ■ 

Proof See Appendix |DJ] □ 

Because w has the same property that Wt + F{w)x — Dw^x G C([0, T], iS(IR)), we 
have the following lemma by repla cing u by w in Lemma |5.7[ 



Lemma 5.11. Let w in (5.15). Given Lq = xq + 7T + in Lemma 
satisfying Assumption 



5.1 



and h''£^ defined by Wt + F{w)x — Dwxx — ^Lq^Lq' 



5.10 



1 



htitr)\\hdt 



i.^o 1 



\\{wt + F{w)x ~ Dwxx)/(t>o{-/Lo)\\l2dt. 



The rest is to compute \\{wt + F{w)x — Dwxx) / 4>a{' / LQ)\\\2dt. Note that 

Wt + F{w)x - DWxx = 7j^[U{x -"fT - Xq) - U{x)] + F{w)x - DWxx- 
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With Lq = 7T + xq + Cw given in Lemma |5.10[ 

\\{wt + F{w), - Dw^,)/M-/Lo)\\hdt 



< 



y2 



< 



2^2 







\\[U{- - - xo) - U]/M-/Lo)\\l2dt 



\[F{w),- Dw,,]/M-/Lo)\\hdt 



{2\\[U{--jT~xo)~U]\\l,+2)dt 



2 1 {2\\Fiw),-Dw,,\\l2+2)dt. 
10 



We note that for |a;o| large, \\U{ — ^T—xq) — 1/11^^2 = Q{\xo\), and \\F{w)x—Dwxx\\'i2 



0(1). Then for Lq = + xq + Cw given in Lemma 5.10 we have the upper bounds 
for J {A) for ja^ol large: 



j{A)<i{w) < e(|xo|/r) + e(r). 



(5.18) 



Now we consider the lower bound. Wc know that (5.13) still holds for |xo| large 
and with Cw given in Lemma 5.10 For |a;o| large, ^[U{x — .tq) — U{x)]dx — Q{\xq\) 
so 



(5.19) 



Consequently, by combining (5.18) and (5.19), we have the bounds in (5.3) and (5.4| 
for \xq\ large: 



J{A) 



>o J 8(|a;o|), for |a;o| large and any fixed T, 

1 9(|a;o|/T), for |xo| large and T small. 



6. Large deviations for discrete conservation laws. Conservation laws can 
only be solved numerically, in general, so we need to consider space-time discretiza- 
tions. For the calculation of small probabilities of large deviations, we may pass 
directly to the calculation of the infimum of the rate function, which we know ana- 
lytically. First we discretize it in space and time and then use a suitable optimization 
method to find the approximate minimizer and the approximate value of the rate func- 
tion. This way of calculating probabilities of large deviations has been carried out 
previously [9, for different stochastically driven partial differential equations. More 
involved methods that use adaptive meshes are discussed in [25^. 

In the next subsection we show briefiy that the rate function for large deviations 
of discrete conservation laws using Euler schemes is, as expected, the corresponding 
discretization of the continuum rate function. 

6.1. Large deviations with Euler schemes. To formulate the discrete prob- 
lem, we discretize the space and time domains with uniform grids, with L = xq < ■ ■ ■ < 
XM = R, MAx = R- L and Q ^ < ■ ■ ■ < tjq ^ T , NAt = T. Let Q'^ denote the 
average of u over the m-th cell [xm-i,Xm] (whose center is a;m-i/2 = i^m-i +a^m)/2) 
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at time nAt, and Q" denote the vector (Q", ■ • • tQm)^- The Euler scheme for the 
conservation law is 



D- 



At 



(Q:;+i-2Q;^ + Q;:._i)+eAW^™, (6.1) 



for m = 2, . . . , M — 1, n = 0, . . . , iV — 1, where Fm±i/2(Q") numerical fluxes con- 
structed by standard finite volume methods such as Godunov or local-Lax-Friedrichs 
(LLF), and possibly with higher order ENO (Essentially No n-Oscillatory) reconstruc- 
tions. The initial conditions (Qm)m=i,...,j\/ are given by (6.4|, the boundary conditions 
{Qi)n=i,...,N and {Q^j)n=i....,N are given by (6.7), and the fluxes are given explicitly 



in ( |6.6[ ), in the next subsection, for Burgers' equation. We let 

W)-(fc2(Q"),...,6A./-l(Q"))^, 

where 



Let {AW^)m=2....,M-i,n=i,...,N be Gaussian random variables with mean and co- 
variance 



E{Aw:^,\Aw::,i) 



At ^ 

'-^mim2 I 



0, 



ni = 7i2, 
otherwise. 



The matrix C — {Cij)ij=2 M-i is symmetric and non- negative definite. For sim- 
plicity, we assume that C is positive definite, and then C ~ for an invertible 
matrix <&. 

By the Markov property, the joint density function /(qo,...,qn) is 



N-l 



Tt = 



where 



/q"+i|Q 


„(g"- 




1 




AtAx 




exp 


2e2 


1 




AtAa; 




cxp 


2e2 



1 



(?"+ - 9") - K?") C 



At 



where = (27re^At/Ax-)*^"-^ detC and the ^^-norm is here 



E 

J=2 



The exact probability that Q — {Q^, . . . , Q^) € A with the given is therefore 

e^) 



ZN 



cxp 



AtAx 
2e2 



JV-l 

E 

n=0 



Ai 



(<z"+^ - g") - 6(9") 



dg^ ■ ■ ■ dq 



N 
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Because the problem is finite-dimensional, the LDP is a form of Laplace's method for 
the asymptotic evaluation of integrals [53] and the rate function for q~ {q^,. . . ,q^) 



IS 

N 



(6.2) 



2 

From this expression it is clear that the rate function of the continuous conser- 
vation law is the limit of (6.2) as At, Ax — >■ 0. However, the LDP of the discrete 
conservation law does not immediately imply that of the continuous case. The limit 
of e and the limit of At, Ax need to be interchangeable. In the language of large 
deviations, the law of the discrete conservation law has to be exponentially equivalent 
to that of the continuous one 0. Without going into the proof of this, we can only 
say that (6.2) is a discretization of the rate function of the continuous conservation 
law. 

6.2. Numerical calculation of rate functions for changes in traveling 
waves. In the numerical simulations we consider the rare event that a traveling wave 
at time becomes a different traveling wave at time T due to the small random 
perturbations. 

We carry out numerical calculations with Burgers' equation as a simple but rep- 
resentative conservation law. For convenience, given a traveling wave Uo with the 
speed 7, we consider Burgers' equation in moving coordinates: 

ut+(^{u--ff^^^{Du,),+eW. (6.3) 



Thus Uq is a traveling wave of (6.3) with speed zero. We are interested in the rare 
event that m(0,x) = U[){x) and u{T,x) = Ut{x) because of eW, where Uq and Ut are 
two different traveling waves. 

We therefore minimize the rate function ( |6.2[ ) to compute the asymptotic proba- 
bility. The initial conditions 

q°, = Uo{xrn-i/2) for all m = 1 , . . . , M (6.4) 

and the terminal conditions 

Qm = UT{xm~i/2) for aU TO = 1 , . . . , M (6.5) 

are simply the values of the functions at the centers of the cells. For m ~ 2, . . . , M —1 
we let numerical fluxes for (u — 7)^/2 at x„i±i/2- We use Godunov's 

method [T7| to construct 

rpn _ /mS;\-i<9<C„ Ul ~ 7)^ C-l < (f. r.^ 

When the final traveling wave solution Ut is the shifted initial profile Uo{x — xq) 
for some xq, then we impose the boundary conditions 

q" = u_ and q'^j = u+ for all n = 1, . . . , iV, (6.7) 

where u± ~ lim^j-^ioo Uq{x). We will discuss the boundary conditions imposed in the 
other cases later on. 
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The covariance matrix of the random coefficients (AWj")i=2,...,M-i is set to be 
equal to Cij = 5ij and thus <i> — /m-2xM-2- In other words, (AWj")i=2,...,A/-i are 
i.i.d. Gaussian random variables. Since the discrete problem is finite dimensional, $ 
is clearly a Hilbert-Schmidt matrix. 

The objective is to minimize the rate function (6.2 1. Because the initial, terminal 



and boundary conditions are easily integrated into the definition of the discrete rate 
function, we can minimize I{q) by unconstrained optimization methods. The BFGS 
quasi-Newton method |19| is our optimization algorithm. 

An important issue in numerical optimization is that any gradient-based method, 
for example the BFGS method that we use, only gives a local optimum unless the 
objective function is convex. In our case, it is not clear that the discrete rate function 
(6.2) is convex or not. However, based on our numerical simulations we note the 



following. 

1. Our numerical results of the optimal shifted profiles coincide with the ana- 



lytical predictions (5.4 1 



2. Instead of using a good initial guess for the minimizer in the numerical op- 
timization, we have also numerically verified that completely random initial 
guesses give essentially the same result. 

3. We have checked numerically to see if the rate function (6.2) is convex. We 



randomly pick two close by test paths to see if the midpoint convexity of ( 6.2 1 



is satisfied. We find that in 99.98% out of 10^ pairs the discrete rate function 
passes the convexity test. It is not known what causes the 0.02% failures. We 



conclude that based on numerical calculations (6.2 ) is essentially convex, if it 



is not fully convex. This explains the observed robustness of the numerical 
optimization. 

6.3. The numerical setup. We use different Uq and Ut to calculate probabil- 
ities of several rare events. Our main interest is the anomalous wave displacement, 
which is theoretically analyzed in the previous sections. The other cases are the wave 
speed change, the transition from a strong shock to a weak shock, and the transition 
from a strong shock to a weak shock. We have not carried out an analysis of the 
last three cases. However, the probabilities of these rare events can be calculated 
numerically and show how unlikely such events are compared to anomalous wave 
displacement. 

In each configuration wc consider the high viscosity case {D = \) and the low 
viscosity case {D = 0.01). T = 1, Ax = 0.2 and At = 0.02 in all simulations and the 
linear interpolation of C/q and Ut is the initial guess in the numerical optimization. 
As we noted before, a random initial guess gives essentially the same result, but we 
use the linear interpolation to speed up the optimization. 

6.4. Anomalous wave displacements. In this case, we let Uq be the traveling 
wave solution for Burgers' equation (6.3), and Ut = Uq{x — Xq) represent a shifted 
traveling wave. This is the discrete version of (5.11, and we find that the numerical 



results are consistent with the analytical result (5.4) 



We first consider Fig 6.1 and Fig 6.2 with D ~ 1. The optimal path is close to the 
linear shift when xo is small and it looks like the linear interpolation when xq is large. 
This also motivates us to choose the test functions v and w for the upper bounds of 



J {A) in Section [O] and [O 

near xq = and is linear for xg large, and is of order 1 /T in T 



Further, Fig 6.2 shows that the optimal / is quadratic 

These observations 



confirms the analysis (5.4) 
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Fig. 6.1. The optimal paths and their values of I in ^6.1^ in the case that Ut{x) = 
Uo{x — xq) for different xq. In each figure, we plot the curves indicating the optimal path at time 
0, At, 2 At, . . . ,T = 1. We can see that for xq small, the optimal path is close to the linear shift 
traveling wave and for xq large, the optimal path looks like the linear interpolation of Uo and Ut- 
These observations are consistent with the test function v and w we use for the upper bounds in 
Section lsTB and 1 5. 31 



We consider next Fig |6.3| and Fig |6.4| with D — 0.01. As the transition regions 
are very narrow and separate very quickly in the low viscosity case, the optimal path 
is nearly the linear interpolation, and therefore the / versus xq plot is almost linear 
except around xq = 0. Moreover, the optimal rate function is still of order 1/T in T, 
which is also seen in the analysis (5.4). 



6.5. Change of wave speeds. In this subsection we consider the rare event 
in which u{Q,x) — Uo{x) and u{T,x) — Ut{x) :— Ui{x) because of eW, where 
Ui{x) is a traveling wave with ui± :— lima;_j.-|-oo Ui{x) which are different from u± = 
limj.^±oo Uq(x), and such that 71 = f(ui+)-f(ui-) ^ different from 7 = -^^"+)~-^^"-) . 
This case does not belong to the class of problems addressed in the previous sections 
of this paper, in which the boundary conditions are fixed. But it is still possible to 
look for the optimal paths going from Uq to Ut and that minimize the rate function 
I. 

The boundary conditions are more delicate to implement in this case. We know 
that Uq and Ut are very close to constants when L and R are far from their transition 
regions. In this case, the LDP implies that the optimal path around the boundaries 
should be very close to the linear interpolations of Uq and Ut in time. Therefore we 
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Wave Shift, D = 1, r = 1, ii_ = 2, u+ = 1 



Wave Shift, D = 1, x„ = 20, u_ = 2, ii+ = 1 





Fig. 6.2. Left: The optimal values of I in \6.2\ ve rsus xq in the case that Ut(x) = Uo{x — xq) 
for xq = 0, 1, 2, . . . , 20 with the same setting in Fig \6.1\ We see that the curve is quadratic for xq 
small and linear for xq large (see ^5^). Right: The reciprocal of the optimal values of I in 
versus T in the case that Ut{x) = Uo{x — xq) for xg = 20 and T = 0.1,0.2, . . . , 1 with the same 
setting in Fig \6. 1\ We see that the curve is linear in T, which means the optimal I is of 0(1/T) 
(see also \5.4\ )■ 



let: 

= (1 - '-^hl + ^-Zf , qM = (1 - + 

It is possible not to set the boundary conditions and to optimize the boundary cells 
as well. Our numerical simulations show, however, that the solution in both cases is 
basically the same except for some oscillations near the boundaries. The oscillations 
come from the inappropriate discretization at the boundaries, and this is a limitation 
of the numerical discretization. We can have a few extra boundary conditions to 
reduce the unwanted oscillations at the boundaries. For example, we can additionally 
set 

q2 = (1 - + ^q^. <1M^, = (1 - + ^-Zfi-i. 



The results are shown in Fig |6.5[ We let 7 = and 71 = 3.5 while we keep 
w_ — u+ ~ wi- — ui+ = 1 to indicate that Uq and Ut have roughly the same transition 
magnitude. We see that the values of the rate function are much larger than that 
of the anomalous wave displacements. This means that it is very unlikely to have 
changes of wave speeds compared to wave displacements. 

6.6. Weak shocks to strong shocks and strong shocks to weak shocks. 

We also consider the case that Uq is a weak (strong) shock while Ut is a strong (weak) 
shock. By a strong shock we mean that the difference between u- and u+ is large. 
This case is also not in the range of our analytical framework, but we can still compute 
the rate function after we impose the suitable boundary conditions. We use the same 
boundary conditions as the ones in the previous subsection: 

q- = {l- ^)g? + ^gf , ql, = {I - ^)q<i, + ^q^„ 
12 = (1 - ^)q'2 + ^1^, IM-I = (1 - + 



From Fig |6.6| and Fig |6.7| we see that the optimal path of weak to strong and the one 
of strong to weak are significantly different, even if the reference strong and weak 
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Fig. 6.3. The optimal paths and their values of I in ^6.1^ in the case that Ut{x) = 
Uo{x — xq) for different xq. In each figure, we plot the curves indicating the optimal path at time 
0, At, 2 At, . . . , T = 1. Different from Fig \6. 1\ as the transition region is very small, we see nearly 
no shifted part but only the linear interpolation. 



shocks are fixed. We note the very large value of the rate function compared to 
anomalous displacement, and how it depends on D. This confirms quantitatively 
the expectation that shock profiles are very stable and they are not easily perturbed 
except for displacements. 

7. Direct numerical simulations with importance sampling. The large 
deviations probabilities calculated in Sections [S] and [6] are only the exponential decay 
rates of the probabilities but not the actual probabilities. In this section we use Monte 
Carlo methods to compute the actual probabilities numerically. 

7.1. Burgers' equation vifith spatially correlated random perturbations. 

We reformulate the discretized problem for Burgers' equation when we have spatially 
correlated random perturbations. Given a traveling wave solution Uo{x — jt) of Burg- 
ers' equation with lima;_>.±oo Uo{x) = u±, we transform to moving coordinates 



ut+(^^iu~jf^^^{Du,).,+eW. (7.1) 



Then Uq{x) is a stationary traveling wave of (7.1 ). The rare event we consider is 



As^{ue £^ such that u(0, •) = Uq, \\u{T, •) - Uo{- - xo)\\l^ < S}. 

Although for a discrete conservation law it is also possible to consider the other cases 
in Section |6j their probabilities are too small to compute by the basic Monte Carlo 
method so we omit them. 
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Wave Shift, D = 0.01, T = 1, u_ = 2, 11+ = 1 



Wave Shift, D = 0.01, lo = 20, ti_ = 2, ii+ = 1 





Fig. 6.4. Left: The optimal values of I in \6.2\ versus xq in the case that Ut(x) = Uo{x — xq) 
for Xq = 0,1,2, ...,20 with the same setting in Fig \ 6. 3\ As what we see in Fig \6.3\ the optimal 
path is the linear interpolation and thus the curves is almost linear except a small perturbation 
around xq = 0. Right: The reciprocal of the optimal values of I in ve rsus T in the case that 

Ut{x) = Uo{x — xo) for T = 0.1, 0.2, . . . , 1 with the same setting in Fig \6.3\ We see that the curve 
is linear in T, which means the optimal I is 0{1/T). 



1=42.4606. D=l 



1=43.0052. _D=0.01 





Fig. 6.5. The optimal paths and their values of I in 6.2 ) in the case that Uo and Ut have 
different speeds. Here we let 7 = 1.5 in Burgers' equation 6.3), and therefore the speed of Uq is 
zero and the speed of Ut is 3.5. The differences between the left and right boundary values are kept 
the same so that Uo and Ut have the same transition magnitude. In each figure, we plot the curves 
indicating the optimal path at time 0, At, 2 At, . . . ,T = 1. 



To compute ¥{u e As) numerically, we discretize the space and time domains 
uniformly as in Section 6.1 L — xq < ■ ■ ■ < xm = R, MAx = R — L and = to < 
■ ■ ■ < tj^ = T, NAt = T. Here denotes the average of u over the m-th cell at 
time nAt, and evolves by the Euler method: 

Q''r;^'=QZr-^iKn+^,-K.-l)+D-^{Ql,^,~2Ql^^ (7.2) 

AIX 2 2 {Ax) 

for TO = 2,..., A/ — 1, n = l,...,iV, where F^j^,^^^{Q'^) are numerical fluxes of 

{u — 7)^/2 constructed by Godunov's method and {AW^)ra=2,...,M-i,n=i,...,N are 
Gaussian random variables with mean and covariance 

'"^ I 0, otherwise. 
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1=524.4324, D=l 1=526.6572, L>=0.01 




-10 -5 5 10 -10 -5 5 10 



X X 

Fig. 6.6. The optimal paths and their values of I in g[ ) in the case that Uo is a weak 
shock while Ut is a strong shock. Here we let Uo has a small difference in the boundary values 
while Ut has a large one. In each figure, we plot the curves indicating the optimal path at time 
0,Ai, 2At, ...,T = 1. 




Fig. 6.7. The optimal paths and their values of I in lj6.S^ in the case that Uo is a strong 
shock while Ut is a weak shock. Here we let Uq has a large difference in the boundary values 
while Ut has a small one. In each figure, we plot the curves indicating the optimal path at time 
0,At,2At,...,T = 1. 



In order to make the problem more realistic, we assume that the variables AW^ are 
spatially correlated: C^jm^ = a'^ ex-p{— j^\xmi — a^msl) ^-nd C = (Cij)fjC2 = 
Finally we impose the initial and boundary conditions: Q'^^ = C/o(a;m- 1/2)7 Qi = u- 
and Q\j = u+. 

7.2. Introduction to importance sampling. To estimate P(m G Ag), we 
may use the basic Monte Carlo method. The Monte Carlo strategy is as follows. 
We generate K independent samples l\W^^^ = {^Wm^'''')m=2....,M-i,n=i,...,N , k = 
1, . . . ,K, of the Gaussian vector {AW^)m=2,....M-i,n=i,....N, which give K indepen- 
dent samples Q^''^ = {Qm^''^)m=i,...,M,n=i,...,N, k = 1,...,K of the random vector 
Q — {Qm)m=i,...,M,n=i,...,N- The basic Monte Carlo estimator is 

P'''' - ^EUAQ''^) (7.3) 

k=l 
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where 

M 

QeAsii and only if Ax ^ [Q^ - Uq{x^_i - xo)]^ < 5^. (7.4) 

m— 1 

In other words, P^^^ is the empirical frequency that QC^) e As. It is an unbiased 
estimator E[P*^'-^] = P(Q e As). By the law of large numbers, it is strongly convergent 
pMC _^ p^Q ^ ^^-j almost surely as K — > oo. Its variance is given by 

Var(P^^^) = lvar(U,(Q)) = -i(P(Q e As) - P{Q e As)'). 

In order to have a meaningful estimation, the standard deviation of the estimator and 
G As) should be of the same order. Namely, the relative error 

Var3(pA^C) If 1 ' 



e As) Vk\V{Q^ As) 

should be of order one (or smaller) . This means that the number K of Monte Carlo 
samples should be at least of the order of the reciprocal of the probability P{Q G As). 
We note that for e small, P(Q G As) decreases exponentially and so K should be 
increased exponentially; the exponential growth of K makes the basic Monte Carlo 
method computationally infeasible. 

The well established way to overcome the difficulty of calculating rare event prob- 
abilities is to use importance sampling. The problem with basic Monte Carlo is that 
for small e there are very few samples in As under the original measure P so the esti- 
mator is inaccurate. In importance sampling we change the original measure so that 
there is a significant fraction of Q'^ in As under this new measure Q, even for small 
s. Since we use the biased measure Q to generate Q^*^-*, it is necessary to weight the 
simulation outputs in order to get an unbiased estimator of P{Q E As). The correct 
weight is the likelihood ratio since we have: 

P(Qe A,) =Ep[U,(Q)] -Eq[U,(Q)^(Q)]. 

Then the importance sampling estimator is 

k—l 

where Q^*^^ are generated under Q. The estimator P^"^ is unbiased Eq[P-^'^] — ¥{Q S 
As) and its variance is 

1 dV 
VarQ(P^«) ^ _VarQ(U,(Q) — (Q)) 

The main issue in importance sampling is how to choose a good Q to have a low 
VarQ[lyi^((3)^(Q)]. In many cases (see for example, [551 HIl 131 HO]), it can be 
shown that the change of measure suggested by the most probable path of the LDP 
is asymptotically optimal as e — >■ 0. However, it is also well-known that in some cases 
(see [12])j the estimator by this strategy is worse than the basic Monte Carlo, and 
may even have infinite variance. However, because the rare event As is convex and 
the discrete rate function / is (numerically tested) essentially convex, the importance 
sampling estimator by this strategy is expected to be asymptotically optimal and we 
will see that it indeed works very well. 
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7.3. Importance sampling based on the most probable path. In this 
subsection we implement the importance sampling by using a biased distribution 
centered on the most probable path obtained in Section [6j From Section 6.1 Q .= 
arginfggyi^ I{Q) is the most probable path as e — > by the large deviation principle. 
We choose h = (h"^)n=i,...,N = {Kn)m=2,....M-i.n=i,....N such that 



D- 



At 

Ax 
At 



(7.6) 



(Ax) 

for m = 2, . . . , M - 1, n = 1, . . . , TV. 
Assume that under the probability 



2(Q™+i-2Q:^, + Q:Ui) + (*n™, 



the vector AW" 



{Aw;; 



in (|7.2| is multivariate Gaussian A/'(£-i$/i", ^C) and Cov(AT^^i , AVF^^J 



if n~ n2. Denoting AW = (ATy")„^i^ 
computed explicitly: 



, the likelihood ratio 



m=2....,Af-l 



can be 



( 1 

(AT4/)=exp --^[|l$-iAW^"| 

\ ri=l 



\^-^ AW" - F-^K 



if;«ii2 



(7.7) 



Note that under Q, (7.2) can be written as 
At 



o"+i = o" 

where AW 



-(F 



1 1 1 , 



-D- 



At 



Ax 

- {AW^)m=2,...,M-l.n=l 



(Ax) 



r(Q™+i-2g™+Q™-i) + ($^")™+eAW^;^, 

(7.8) 

is zero-mean, Gaussian with the spatial 



covariance -^C, and white in time. Then (7.7 1 can be written as 

^(AW^) = exp l^-i J2 ll'&^'AVK" + e-'hm - \\'f-'AW"\\l^ 



(7.9) 



In summary, the importance sampling Monte Carlo strategy is implemented as follows: 

1. Compute the optimal path Q = arg infgg^^ /(Q) and its residual h from 

2. Sample K independent AW^'^'^ with the zero-mean, Gaussian distribution 
with the covarian ce in space and white in time. Compute the corre- 
sponding QC^') by (frSl. 



3. The importance sampling estimator is 

K 

K 



fc=i 



where ^{AW) is defined in (7.9). 

7.4. Simulations with importance sampUng. We consider three estimators: 
the basic Monte estimator P^ ^*^ and two importance sampling estimators: Pn^ and 
P^^, where P^'^ uses h in ( |7.6| with Q — arginfgg^^ HQm) while Pg^ uses h in ( |7.6| 
with Q = arg inf gg^^ I{Q)- The parameters for the simulations are listed in Table 

rm 
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Ax = 0.5 


At = 0.05 


L = -15 


R^20 


r = 1 


u_ = 2 


= 1 


7 = 1.5 


D = 1 


xo = 5 


S = VO.5 


cr = 1 


lc^5 


iiT = 10* 



Table 7.1 

T/ie parameters for the Monte Carlo simulations. 



inf„E,4„/(ii) = 0.093421 int„gA, = 0.017911 




Fig. 7.1. T/ie optimal paths for vaiq^Ao I{Q) andiniQ^A^ I(_Q)- With the spatially correlated 
noise, inlgg^^ I{Q) much lower than the one with the spatially white noise (see Fig \6. 1^ . 

As we mention before, we only test the wave displacement with = 5 and D — 1 
because the probabilities of the other cases in Section [g] are too small for p^'^'^ to have 
meaningful samples. 

First we find the optimal paths for infgg^g I{Q) ^^nd infgg^^ I{Q)- As before, 
infggyip I{Q) can be modeled as an unconstrained optimization problem and we solve 
it by the BFGS method while we solve infgg^^, I{Q) by sequential quadratic program- 
ming (SQP) nn]. 

From Fig |7.1| we note that inlgg^^, I{Q) is much smaller than the corresponding 
one with spatially white noise. This is because with the correlated noise, it is easier 
to have simultaneous increments. In addition, infgg^^ I{Q) is significantly different 
from infgg^Q I{Q) because i5 = v^OTS is not very small, we will see that this difference 
significantly affects the performances of P^^ and Pg^ . 

Once arginfgg^,-, I{Q) and arginfgg^^, I{Q) £^re obtained, we can construct Pq^ 
and P^^ ■ We estimate P(Q G As) by these three estimators for 100 different e taken 
uniformly in [0.01,0.2]. For each e and estimator, we use K = 10* samples. The 
(numerical) 99% confidence intervals are defined by [Mean if — 2.6 Std^, Mean + 
2.6Stdif] where the (numerical) mean and standard deviation are 

Mean,, = 1 f^pC^), Std^ = 1 Y.{p^'^f Mean^^, 

k=l k=l 

with pC^) = 1a,{Q^^^) for P*''^ and pC^'' = dF / dq{AW^^^)lAAQ^^^) for P^^ . We 
find that P/'^ has the best performance and Pq^ also has the good performance for 
0.1 < £ < 0.2. For e < 0.1, because arginfgg^Q I{Q) is not the optimal path, P^^ is 
even worse than P^^'^ due to the inappropriate change of measure. For e < 0.1, be- 
cause arginfgg^j I{Q) is the optimal path, Pg^ dramatically outperforms P'^^^ . This 
shows that large-deviations-driven importance sampling strategies can be efficient to 
estimate rare event probabilities in the context of perturbed scalar conservation laws. 
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Fig. 7.3. Sample paths under these three probability measures with e = 0.1. 



We plot the estimated probabilities and the relative error, the ratio of the (nu- 
merical) standard deviation to the (numerical) mean, in the log scale. Note that in 
the extreme case, the estimated probability is dominated by the value pC^o) of one 
reahzation over K (pC'") = 1 for P^^^- ^nd pC^o) = dP/dQ(AW^(*=o)) for P^-^), and the 
numerical variance is approximately {p^^^^)"^ / K . Therefore the r elati ve error is about 



7.4 



saturate at 100, 



'K = 100. This is why the curves of the relative errors in Fig 
and it also tells us that when the relative error reaches \/K, the estimator is in the 
extreme case and so is not reliable. 

8. Conclusion and open problems. We have analyzed here the small proba- 
bilities of anomalous shock profile displacements due to random perturbations using 
the theory of large deviations. We have obtained analytically upper and lower bounds 
for the exponential rate of decay of these probabilities and we have verified the ac- 
curacy of these bounds with numerical simulations. We have also used Monte Carlo 
simulations with importance sampling based on the analytically known rate function, 
which is efficient and gives very good results for rare event probabilities. 

Acknowledgment. This work is partly supported by the Department of En- 
ergy [National Nuclear Security Administration] under Award Number NA28614, and 
partly by AFOSR grant FA9550-11-1-0266. 

Appendix A. Proof of Lemma |2.2[ Taking a Fourier transform in x we have 

Jo 

where W(t, £,) — J W{t, x)e~^^^dx is a complex Gaussian process with mean zero and 
covariance 



E[#(t,e)M^(i',e)] =tht'C{^,0, 
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Fig. 7.4. The semilog plots of the estimated probabilities (left) and the relative errors (right). 

with C{^,C) = / e-'^''+'^''''Cix,x')dxdx'. We find that 

On the one hand, the fact that $ is Hilbert-Schmidt from into implies that 

(7(e,orfe = 2^ll<f(-,-)lli^(ExE) 



is finite. On the other hand we have (1 — e '^)/s < 1 uniformly with respect to 
s G (0, oo), so we get that for any t e [0, T] 

'(l+aE[|^(t,OI']rfe<27r(r+^)|l<i>(.,.)lli2(RxK), 
which is equivalent by Parseval relation to 

nwzmu < (T+^)ii<i>(-,-)iii.(RxB), 

which gives that Z{t) e L^([0, T], iJ^(M)) almost surely. Similarly we find 
E[|Z(i, - Z{t',0\'] = [1 - e-^«'l*-*'l] [2 - e-2^«'*^*' + e-^«'(*+*')] 
<[l-e-^«^l*-'l]^<|i-t'|C(C,0, 

which gives 

E[\\Z{t) - Z{t')\\l2] = JJ E[\Zit,x) - Z{t',x)f\Z{t,x') - Z{t',x')\'']dxdx' 
< JJ E[\Z{t,x) - Z{t' ,x)\^]^^\[\Z{t,x') - Z{t' ,x')\''Y^^dxdx' 
= jj E[\Z{t,x) - Z{t' ,x)\^]E[\Z{t,x') - Z{t' ,x')\'^]dxdx' 
= 3( / E[|Z(t,x)-Z(i',x)|2]dx)'<3(2^)2|t-if 



26 



where we have used the fact that Z(t, x) — Z{t' , x) is a Gaussian random variable and 
Parseval equahty. By Kolmogorov's continuity criterium for Hilbert valued stochastic 
processes [4, Theorem 3.3] we find that that Z(t) € C{[0,T], L'^{R)) almost surely. 

If $ is Hilbert-Schmidt from into , then 



eCitOd^^27r\\d,^;-)\\l. 



is finite, and we can repeat the same arguments to show the final result. 
Appendix B. Proofs in Section [3[ 



B.l. Proof of Propsition 



3.2 



2.2 



Z{t) = /q S{t - s)dWis) is 



From Lemma 

in C{[0,T], H-^{M.)) almost surely. We want to prove the existence and uniqueness in 
£^ of the solution to the equation 



u(t)^S{t)U- / S{t- s)[F{u{s)%ds + eZ{t), 
Jo 



(B.l) 



where is a function with max{||i^'||i=o , ||_F"||i=o} < Cp < oo. 

To prove the existence we use the Picard iteration scheme: we define = U 

and 



,Tl+l 



{t) = S{t)U- I S{t- s)[F{u'\s))]^ds + eZ{t). 





Therefore, 



S{t-s)[F'{u-{s)){u-{s)-U,)]ds 



(B.2) 



S{t - s)[F'(u"(s))C/:r]ds + eZ{t). 



It is easy to see that — U E C{[Q,T], H^{M.)) for all n, and we first show that 

supjg[o,T] ^ ^ll_ffi(R) a-re uniformly bounded in n. 

Lemma B.l. Let a„(i) = — U\\hi. Then there exists a constant Cu such 

that Qnit) < Cu for all n and t G [0,r]. As a consequence, we can choose sufficiently 
large Cu such that ||u"(t)||/^2 < C„ for all n and t € [0,T]. 

Proof Note that sup(g[o.T] W'^WU — U\\h^ £^nd supjg[o,T] are bounded. 

In addition, by Chap. 15, Sec. 1], \\S{t ~ s)\\c(l^ ^h^) = 0{{t - s)-^^^), then 



5(t-s)[F'(w"(s))[/^]ds 



< / \\Sit ~ s)\\ciL-,Hi)\\F'{u''{s))U,\\L2ds 

ffi Jo 

Jo 
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are uniformly bounded in n, and 



S{t-s)[F'{u-{s)){u:{s)^U,)]d^ 



< / \\S{t-s)\\ciL-^m)\\F'{u-{s)){u:is)-U.,)U2ds 



<Cp / \\S{t- s)\\c(L-^m)\\ul{s) -U^^L-ds 



<Cf / \\S{t^s)\\ciL-^m)\\u'\s)-U\\Hids 



< Cf 



ds 



■n 1/P 






[i 



\\u-{s)-ur„,ds 



1/9 



By letting l/p + l/q = 1 with 1 < p < 2 and 5 > 2, we can have the following recursive 
inequality with a sufficiently large C: 



an+i{t) < 



C C 
2 + 2 



al{s)ds 



1/9 



By the convexity of a; 1— >■ x', 



a^+iW< f + f 



al{s)ds 



uo 



< ^ + ^ I aU.s)ds. 



(B.3) 



By noting that ao{t) = and (B.3), it easy to see that a'^{t) are uniformly bounded 
in n and t S [0, T] and so are a„(t). □ 

To prove the convergence of u"{t) — J7 in C([0, T], H^{^)), it suffices to prove that 

Er=oSUPte[0,T] IK^Ht) -U"{t)\\m < OO- 
LEMMA B.2. Letbnit) = Then Y,'^^QSupte[o,T]bnit) < 00. 
As a consequence, u^it) — U converges to u[t) — U in C([0, T], iJ^(IR)) as n ^ 00 and 
u{t) solves iB.l ). 



Proof. By noting (B.l I, Lemma B.l and || • ||l°° < Cs\\ ■ \\m (the Sobolev embed- 
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ding), we have 



Jo 



+ / \\S{t^s)\\ciL-^m)\\u:-\s)[F'{u-is))~F'iu-\sm\L^ds 



\\S{t~s)\\c(L^^m)\\u'\s)-u^-\s)\\Hids 



+ / \\S{t-s)\\c,^i^.,H^)\K-\s)U4F'iu-{s))-F'iu--\smL^ds 
Jo 

<Cf f \\Sit- s)\\ciL-.Hi)\\u"{s) - u''-\s)\\Hids 
Jo 



\W'{s)~u^-\s)\\%,ds 



< Cf{1 + CsCu) 



- 3)\\c{L\H^)\\u"is) - 








1/p 




\\Sit- s)\\l^^2^H^)ds 




[/ 

Jo 
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where l/p + 1/q = 1 with 1 < p < 2 and q > 2. Then there exists a constant C such 
that 

bl^.it) < C f bUs)ds. 
Jo 

Then bl{t) < bl{t)C''T'^/n\ so b,,{t) < boit){C"'T'' /niy/^ and it is easy to see that 
Er=oSuPte[o,T] bnit) < oo. □ 

Finally we show that (B.l) has a unique solution in 

Lemma B. 3. If u,v G solve [B.l], t/ien supjgjo t] ll""!^) ~ ^(Oll-ff^ = 0- 
Proof. Let m, u G £^ solve (B.l). Using the same calculations in the proof of 
Lemma [B.2[ we get 



\\u{t)-v{t)\\H^ <Cf{1 + CsCu) / \\S{t - s)\\c(L-^.^H^)\\u{s) - v{s)\\Hids. 

Jo 

By noting that ||u(0) — i;(0)||/fi = 0, supf^jg yj \\u[i) — v{t)\\Hi = by Gronwall's 
inequality. □ 



B.2. Proof of Proposition [373| The LDP can be obtained following the strat- 
egy of [m [ni . The first step of the proof uses the LDP for the laws of the stochastic 
convolution eZ on the space C([0, T], iJ^(IR)), where 

Z{t) = f S{t-s)dW{s). 
Jo 

The laws of eZ are Gaussian measures and the LDP with a good rate function is 
a consequence of the general result on LDP for centered Gaussian measures on real 
Banach spaces [HI E] ■ 

The second step is to prove that the mapping X{t) u{t) — U is continuous from 
C([0,T],iJi(M)) into itself, where 



{t) = S{t)U- / S{t- s)(F{u{s)))^ds+ / S{t-s)X{s)ds 



(B.4) 
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Then the LDP for is obtained by the contraction principle [7]. 

To prove the continuity of the mapping, given two pairs {X, u) and (F, v) satisfying 
(|R4l), we show that v - u ^ as Y ^ X in C([0, T], ffi(IR)). By using essentiaUy 



the same calculations in Lemma B.2 



\\v{t)^u{t)\\Hi<CF{l + CsCu) I \\S{t-s)\\ciL^-^m-^\\v{s)-u{s)\\H^ds 

Jo 

+ f \\S{t - s)\\c(H\H^)\\Y{s) ~ Xis)\\Hids, 
Jo 

where C„ = sup^gjo j.] \\ux{t)\\L^. Noting that \\S{t - s)\\c{m,m) = 0{1) and \\S{t - 
s)\\c{L^,H^) = 0{{t — s)^^/^), and by Gronwall's inequality, there exists a constant C 
such that for te [0, T] 

Ht) - uit)\\Hi < Ce^fo(^~^r'''^^ r ||y(,) „ X{s)\\H^ds. (B.5) 

Jo 

Because {T - s)-'^l^ds < oo, v - u ^ Q &sY ^ X \n C{[0,T],H^{M.)). 

B.3. Proof of Proposition [3]4[ For each n, we let /i" e ^^([O, T], L2(M)) such 
that u = and 

\ r \\h-it,-)\\ut--<i{u)<i r \\h-it,-)\\ut. 

2 Jo n 2 Jo 

Because u ^ U{x — -ft), I{u) > 0. We let m" be the mild solution of 

< + (FK)), = (i^O, + (1 - 
u"(0,a;) = 

Then u" - [/ e C([0,T],iJi(M)) and 

1 1 1 



2 Jo 2 Jo 

We show that w" - u ^ in C([0, T], iJi(M)). Let a" = (1 - {nl{u)y^y/^ and 
replace (X, u) and {¥, v) by (u, and (u", a"$ft.") respectively in (B.5 1. We have 

/■T 



\\u^{t) - u{t)\\Hi < Ce^'Ioi^-^r'^^'i^ / (1 _ 



Then - in C([0, T], ^^(R)) as a" ^ 1. 



Appendix C. Proof of Proposition |4.1[ 

We use the following technical result: If / e C([0, T],H^{R)), then F{t) = S{t- 
s)dxf{s)ds is such that / F[t,x)dx = for any t e [0,T]. Indeed, since /(s) e i/^ 
and iJ^ C we have by integrating by parts: 

^^(i,a;) = -y^ -^===J dy{e-^)f{s,y)dyds, 
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so that, for any a < b 

rb 



F{t,x)dx= / , ^ / (e ^('-) -e ^i^-^) ) f{s,y)dyds. 

Jo y'2n(t — s) J 

Let S > 0. There exists M such that J^_j^j j^^^^ f{s,yydy < for any s G [0, T]. 
We can split the integral in y into two pieces. The integral over [— M, M]'^ can be 
bounded by Cauchy-Schwarz inequality so we get 

^'5+/ — ===== / (e ^('-)+e ^(^-^^)f{s,y)dyds, 

and therefore, by Lebesgue's theorem 



lim sup I / F{t,x)dx\ < -ji=S. 

— oo,6— >+C30 Ja V ^ 



Since 5 is arbitrary we get the desired technical result. 

We then find by integrating ( B.2[ ) in x that the center of is given by (4.2 ). With 
the condition that C is in L^{M. x M), the last term is proportional to a Brownian 
motion. 

Appendix D. Proofs in Section [5} 

D.l. Proof of Lemma |5.4[ The values of cr, Lq and Ic do not affect the result, 
so in the proof we set a = Lq = Ic = I and <i>^^ = $ without loss of generality. 
^{x,x') is Hilbert-Schmidt from to H'' if and only if dl^{x,x') € L'^{R x R) for 
< j < k. Taking the two-dimensional Fourier transform on d^^{x, x'), we have 

T^,i'{di<i>{x, x')} - T^{T^,{dmx,x')}} - Tadi[M^)T^'{Mx ~ ^')}]} 

= j-aa^[0o(x)e"«'0i(-e')}]} = m'M-emM^V''^'} 

By a simple calculation it is easy to see that (iO"''/'i(~^')</'o(^ + C') G x M) if (f>Q 

and are both in H^{M.). 



D.2. Proof of Lemma 5.5 , It suffices to show the case that (5„ = 1 /n. For each 
n, let u„ e ^i/„, such that J'(Ai/„) < /(u") < + 1/n; {/(u")} are bounded 

from above by J{A\) + 1 < oo. Because / is a good rate function and compactness 
is equivalent to sequentially compactness in f ^, has a convergent subsequence 

{m"' } whose limit u* is in A. As / is lower semicontinuous, then 

J {A) > lim J^(Ai/„) = lim/(u"') = liminf > /(u*) > J (A). 



D.3. Proof of Lemma 5.6, Because (j)o{x) = 1 on x G (—1, 1) by Assumption 

\m\-/Lo) - IPA- - {t/T){jT + xoMh 

-La I'oo 

+ / {[00 i(a;/Lo) - l]U,{x - {t/T){-fT + x^Wdx 

yo J Lq 
-Lo-ipilT+xo) 

{[0^i(x/Lo + 7^(7^ + xo)/Lo) - l]U.,{x)Ydx. 

Lo-^i-yT+xo) ^ 
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Let Lq = 7T + xq + Cv with C„ > 1. Because 0g ^(x) > 1 is monotonically increasing 
for a; > and Cy > 1, 

[<Po\x/Lo + |(7r + xo)/Lo) - If < [(f>o\x + 1) - l]^ a; > 0. 

Similarly, because 4'q^{x) > 1 is monotonically decreasing for x < 1, 

[(t>o\x/Lo + |(7T + xo)/Lo) - 1]2 < [(l)-\x) - x<0. 

In addition, < Lq — ^{jT + xq) and —C^ > —Lq — y (7T + xq). Therefore, 



m^\-/Lo)-l]U,{--{t/T){jT + xom 



2 

L2 



{[00 - l]C/,(x)}2dx + / {[<f>^'ix + 1) - l]C/,(a;)}2dx. 



< 



By noting that i/iq ^ (x) has at most polynomial growth, there exists a uniform Cy > 1 
for ( |5^ . 

D.4. Proof of Lemma |5.7[ It is easy to check from the properties of the trav- 
eling wave U that Vt + F{v)^ - Dv^^ € C([0, T], By and ([5J]), we have 

vt+ F{v)a;- Dva^x = <J<Po[^) ^01 (^) * /^Lq a;) . 
Taking the Fourier transform on the both sides, 

^L„(*'0 = ^'\MlcOr'J'di^'t + F{v)^ - Dv^^)/M-/Lo)}- 

Because Vt + F{v)x — Dvxx € SiM) and l/0o('/io) G C°° has at most polynomial 

growth, F^{{vt+F{v)x - Dvxx)/M-/Lq)} is well-defined and also in C([0, T], 

In addition, l/\4>i{lcO \ ^-l^o has at most polynomial growth and then indeed h'^£ (t, ^) G 

L2(M). 

Because T^{ivt + F{v)x - Dvxx)/M-/Lo)} is in C([0,r],5(M)) and l/|<^i(C)| has 
at most polynomial growth, 

^ ^ ■ li{0)\-^\F^{{vt + F{v)x - Dvxx)/M-/Lo)}\^d^ 



\\hl{t,-)\\h 



27r(T2 
1 



Wivt + Fiv)x - Dvxx)/M-/La)\\ 



L2- 



as Ic —7- 0, uniformly in i e [0,T]. Then we have (5.10) 
D.5. 

\\Ht,-)\ 



D.5. Proof of Proposition pTs) . For any pair {u,h) satisfying (3.1), 



\h = sup 



if J) 



Because <i>^£^h £ C{[0,T], H^(R)), we have u-U e C([0, T], iJi(M)). Therefore 
VLm.x^±oo u{t, x) — u± and lima;_j.±oo "^(i, a;) = for all t € [0, T]. Thus, 



{ut + F{u)x ~ Duxx, l)^dt 



ut{t,x)dx + F{u+) - F{u-)] dt 
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By the Rankine-Hugoniot condition (2.3) 



u+ — u_) j dt 



[ut(t, x) — —U{x — jt)]dx ] dt 



= \mLof^\\l'l (yjJ[u{t,x)-U{x-^t)]dxj dt. 



By the Schwarz inequality and noting that u G A, 



— / [u{t, x) - U{x - "ft)]dxdt 



T-m'flflW-J 



u{T, x) — U{x — jT)]dx 
lT-'m'£yi\\-J {^j [U{x - 7T - xo) - U{x - ^T)]dx 
\T-^m%fl\\-J [U{x - xo) - U{x)]dx'j . 



D.6. Proof of Lemma 5.9, We first compute ($^^)-^. For any test functions 
/ and g, 



x , f 1 , ,x — x' , 



Lq J I. 



= I f{x') / cjM^)l^<l^i{'^)9{^)dxdx' = {f,{'^'lfg). 



Thus {<^%fg{x) = a[</>o(^)g(x)]*^0i(-|^) and {'^\fl{x) = a0o(^)* t^^iC-^)- 
Then 



1 



D.7. Proof of Lemma 5.10, The proof for F(w)x — Dwxx is similar to the 



proof of Lemma 5.6 so we skip it. For U{x — — xq) — U{x), because (poix) = 1 on 
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e (—1,1) by Assumption 5.1 



m^\-/Lo)^l][U{--jT-xo)^U]\\l. 

+ / {[cPo\x/Lo) - l][U{x ~ ~ xo) - U{x)]fdx 



— OO 

Lo 



< / {[<f>^\x/Lo)~l][u_~U{x)]ydx 



OO 



+ / {[(j)o\x/Lo)-l][U{x-jT-xo)-u+]ydx 

JLo 

-Lo 

{[(l>^\x/Lo)-l][u^-U{x)]fdx 

{[q^o'i^/Lo + {iT + Xo)/Lo) - l][Uix) - u+\fdx 

I Lq-'^T-Xq 

{[<PoHx) - - Uix)]Ydx + / {[0o-i(x + 1) - l][Uix) - u+]ydx. 

-OO JC^ 

The last inequality holds because (j)Q^{x) > 1 is increasing for a; > 0, decreasing for 
a; < and Lq > 1. Then we can find a uniform C^, > 1. 
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